Epidemic processes with immunization 
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We study a model of directed percolation (DP) with immunization, i.e. with different probabilities 
for the first infection and subsequent infections. The immunization effect leads to an additional non- 
Markovian term in the corresponding field theoretical action. We consider immunization as a small 
perturbation around the DP fixed point in d < 6, where the non-Markovian term is relevant. The 
immunization causes the system to be driven away from the neighbourhood of the DP critical point. 
In order to investigate the dynamical critical behaviour of the model, we consider the limits of low 
and high first infection rate, while the second infection rate remains constant at the DP critical 
value. Scaling arguments are applied to obtain an expression for the survival probability in both 
limits. The corresponding exponents are written in terms of the critical exponents for ordinary DP 
and DP with a wall. We find that the survival probability does not obey a power law behaviour, 
decaying instead as a stretched exponential in the low first infection probability limit and to a 
constant in the high first infection probability limit. The theoretical predictions are confirmed by 
optimized numerical simulations in 1 + 1 dimensions. 

PACS numbers: 05.70.Ln, 64.60.Ak, 64.60.Ht 



I. INTRODUCTION 



Epidemic processes can be described as the spread and 
decay of a non-conserved agent, an example of which is 
an infectious disease 0. The agent is not allowed to 
appear sponstaneously but it can multiply itself by in- 
fecting neighbouring individuals, or decay at a constant 
rate. Depending on the balance between these two pro- 
cesses, the infection may either die out or spread over 
the entire population. The two regimes of survival and 
extinction of the epidemic are typically separated by a 
continuous non-equilibrium phase transition. When the 
decay process dominates, the epidemic dies out at large 
times and the system gets trapped in an absorbing state 
from which it cannot escape. 

Continuous phase transitions into absorbingstates are 
associated with certain universality classes 0, Q- For 
epidemic processes, a well studied case is the univer- 
sality class of directed percolation (DP), ft is believed 
that two-state spreading processes with short-range in- 
teractions generically belong to the DP class, provided 
that quenched randomness, unconventional symmetries 
and large scales due to memory effects are absent 0, . 
Examples of physical systems whose critical behaviour 
is described by DP include heterogeneous catalysis @, 
chemical reactions @, HI , interface depinning 0, Fnl . the 
onset of spatio-temporal chaos , flowing sand |12| and 
self-organized criticality 01 ■ 

The epidemic process in which the susceptibility to in- 
fection is independent of previous infections is described 
by DP. However, for a more realistic description, wc 
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should consider an immunization effect |l4j . Immuniza- 
tion can be added to the DP model by changing the sus- 
ceptibility after the first infection [lol fill . Il7| . A min- 
imalistic model that captures this feature is one that is 
controlled by two independent parameters: a probability 
of first infection and another probability for all subse- 
quent reinfections. The fact that the local susceptibility 
depends on whether a site has been infected in the past or 
not leads to a non-Markovian epidemic process, in which 
the time evolution depends on the entire history. This 
non-Markovian feature changes the universality class of 
the epidemic spreading. 

The phase diagram of an epidemic process with im- 
munization (Fig. ^) was studied in Ref. PHH. If the 
probabilities for first infections and reinfections are equal, 
the model corresponds to ordinary directed percolation. 
However if the susceptibility changes to zero after the first 
infection there is perfect immunization, and the model re- 
duces to the General Epidemic Process (GEP) 0. GEP 
belongs to the ordinary percolation universality class [2^ • 
The critical points of the GEP and DP are connected 
by a curved phase transition line separating a phase in 
which the spreading process always dies out, from an- 
other phase of annular growth, where an active front 
may propagate into regions of non-immune sites, leav- 
ing a bulk of immune sites behind. As shown in Ref. 0] 
the critical behavior along this line (except for the upper 
terminal point) belongs to the same universality class as 
GEP. Using field theoretic renormalization group tech- 
niques, the critical exponents were calculated along this 
line [!l E3- The main result in Ref. [3 is that the 
compact growth/no growth phase transition line is at the 
critical value of the reinfection rate and independent of 
the first infection rate. Above this horizontal transition 
line in Fig. ^ the model exhibits compact growth and 
approaches the stationary state of supercritical DP. This 
is because, in the active phase of an epidemic process 
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diction is confirmed by optimized numerical simulations 
in 1+1 dimensions. The high first infection probability 
limit is studied in Sec. lIVI giving similar results. We com- 
plete the section with a theoretical approximation for the 
spreading behaviour, and with numerical calculations to 
corroborate the theoretical claims. Finally, a discussion 
of these results together with a possible connection with 
multiple absorbing state models is the subject of the con- 
clusions in Sec. 



II. FIELD THEORETICAL ANALYSIS 



The model 



FIG. 1: Phase diagram of directed bond percolation with im- 
munization in 2+1 dimensions. Along the curve phase tran- 
sition line the universality class corresponds to the GEP. The 
horizontal line separates the no growth-annular growth region 
from the compact growth behaviour. The point where both 
phase transitions lines meet correspond to the universality 
class of DP. 



with immunization, each site will be visited at least once 
after a sufficiently long time so that the dynamics in the 
stationary active state involves only reinfections. On the 
horizontal phase transition line itself all reinfection pro- 
cesses are critical DP while the probability of first infec- 
tions may be sub- or supercritical. By varying the first 
infection rate, we can impede or facilitate the spreading 
into non-immune regions. In Ref. [Tsl ] a numerical analy- 
sis of the scaling behaviour along this horizontal line gives 
no hint of power law behaviour, and it is also suggested 
that this result could be applied to models with multiple 
absorbing states [IliJlHEHmiHlillil. 

The aim of the present work is to investigate in further 
detail the dynamical critical behavior along the horizon- 
tal phase transition line and in the vicinity of the DP 
critical point, and to give a theoretical explanation of 
the absence of power law scaling along this line. 

In order to describe the effect of immunization in an 
epidemic process, in Sec. [H] we study a field theoretic for- 
mulation of directed percolation with immunization and 
we show that the non-Markovian term which contains the 
immunization effects, is relevant under renormalization 
group analysis. As a consequence of this result we argue 
that the asymptotic spreading properties along the hor- 
izontal transition line should be determined by the lim- 
its of very low and very high first infection probability. 
Therefore, in Sec. lIIII we present a study of the very small 
infection probability limit and we develop a quasi-static 
approximation to obtain the scaling behaviour of the sur- 
vival probability of the epidemic. It turns out that this 
does not follow a power law, but instead decays asymp- 
totically as a stretched exponential. This theoretical pre- 



In this section we develop an alternative derivation of 
the action for the DP model with immunization that was 
proposed in [TH fill lri| . The microscopic rules for DP in 
d + 1 dimensions are rather simple. An infected site at 
time t can infect its nearest neighbours at time t + 1 with 
a probability po. There is a critical threshold p c such 
that for po < Pc the epidemic process always dies out, 
that is, it reaches the absorbing state. For po > p c there 
is a finite probability that the epidemics survives. At the 
critical point po — p c the system scales anisotropically in 
time and space. The upper critical dimension is d c = 4, 
below which the fluctuation effects become important. 
The field theoretic action of DP E9,|33 reads as follows: 



5bp = J dtd d x[4>{d t -DV 2 +r)(t)^u 1 4)4) 2 -u 2 4) 2 (f\. (1) 

Here, <j> is the local activity, <j) is the response field and 
r oc p c — po is the mass parameter which measures the 
distance from criticality. This action can also be written 
as a Langevin-type equation for the local activity, 

(d t - DV 2 + r)tf> + X -u<\> 2 + £(x, t) = , (2) 

where u is the symmetrized coupling constant after 
rescaling the fields according to the DP time reversal 
symmetry. The noise £ (x, t) is Gaussian, and satisfies 
(£(x, t)> = and (£(x, t)£(x', t')) = u5 d (x - x')6{t - t'). 
For a systematic analysis of the immunization around the 
DP fixed point it is often more convenient to make use 
of the description in terms of the action. 

To add immunization to the model, the susceptibility 
after the first infection is decreased by an amount A > 

0. The probability of infection, p, depends locally on 
position and time, p = p(x, t). The microscopic rules are 
modified as follows: a healthy site can first be infected 
with infection probability pq. A site which is infected 
at time t will become immune at the next time step t + 

1. Meanwhile any immune site can be re- infected at the 
lower re-infection probability po — A. Thus, the state of 
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the system at time t depends not just on the configuration 
of infected sites at time t — 1, but on the entire previous 
history. 

In order to take into account the effects of immuniza- 
tion, we will modify the DP action. Let us consider a dis- 
crete d + 1 dimensional lattice and assume that the field 
0(x, t) is almost constant between times t and t + At, 
where At = 1 is the time step unit on the lattice. We 
then subdivide it into N intervals, At/N. First of all, we 
want to find an expression for the probability that there 
is an active site in the temporal sub-interval At/N. To 
do so in terms of the field </>(x, t), we need to take into ac- 
count the fact that <f> has units of length~ d / 2 . Then this 
probability can be written as w<p(x, t)At/N, where we 
introduce a parameter w to ensure that this expression 
is dimensionless. 

In this way, 1 — wcf>(x, t) 4£ is the probability that the 
site is not active in the interval At/N. The probability 
of not finding activity between times t and t + At, can 
then be expressed as: 



N 



At. 



Hil-w^t + i/N)— ) 



i=0 



(l-^(x,t)^) w 



N- 



Finally the action for DP with immunization can be 
written as a modification of the directed percolation ac- 
tion as follows: 

S = 5dp + 

A / dtd d x4>(x, i)^(x, t) [l - e~ w it dt >( x >*')] . (7) 

According to the field theory, w is a finite coupling 
constant. Then in the numerical calculations on a dis- 
crete d + 1 dimensional lattice, w can be considered as a 
parameter such that the subsequent infection rate is, 



Pa + (P - Pa) [1 - exp(-wn(x, t))} 



(8) 



-wc/>(x,t)At 



(3) 



Here, the function n counts all the past activity at a site 
x until time t. Since the exponential function vanishes at 
sufficient large times, the effective subsequent infection 
rate is equal to p. 

In the numerical calculations carried out in the present 
work, we assume that the susceptibility for spreading 
only changes at the first infection process and remains 
constant thereafter. This assumption is equivalent to 
take an infinite value of w, but does not make any rel- 
evant change in the final results of the theory as it is 
argued in Ref. [Tsl ]. 



However, between times t = and t > 1, we can no 
longer assume that the field 4> is constant. Therefore, the 
probability of a site x has never been infected by time t 
turns out to be, 



-tu0(x,t')At 



(4) 



When the continuous limit is taken on the discrete lat- 
tice, the probability of having been infected at least once 
in the past becomes 1 — exp(— w J <p(x,t')dt'). This ex- 
pression is the probability for a site to be immune at time 
t. Since A is the parameter related to the immunization 
in the system, we can use the above result to write an 
expression for p(x, t) , 



p(x, t) = po - A 



1 



-w J* <Xx, t')dt' 



(5) 



We assume that the mass parameter r can be written as 
r cx p c — p(x, t). Thus, using Eq.JS] the addition of immu- 
nization can be reflected in the action as a modification 
in the original DP mass parameter by the following sub- 
stitution: 



r + X 



-w J" * 0(x, t')dt' 



(6) 



B. Renormalization group analysis 

Expanding in Eq.[7|the exponential function as a power 
series, the action reads 

S = S D p+ / dtd d x4>(x 7 t)(f>(x 7 t) x 



~ aw 



dt/<f>(x, t') 



(9) 



where we abbreviate A^™) = Xw n . Let us now study the 
stability of the DP fixed point with respect to a small 
perturbation in A*™) by dimensional analysis. Introduc- 
ing a length scale the mean-field fixed point of DP is 
characterized by the dimensions [00] = k d , [D] = uik~ 2 , 
and [u/D] = k^/ 2 so that the upper critical dimension 
is d c = 4. When immunization is considered as a pertur- 
bation, the effective expansion parameter in the expres- 
sion of the two point vertex function is u^X^"'- Thus, by 
dimensional analysis we find 



i D 2n+l- 



k 2+n(4-d) 



(10) 



Consequently the upper critical dimension for each cou- 
pling constant A'-™' is d c = — + 4, where n — 1, . . . , oo. 
This means that for dimensions d > 6, A^™) is irrelevant 
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FIG. 2: One-loop diagrams for the three point vertex function T^^. The dashed lines in the external dashed-continue-leg 
represent the shift in time indicated by the Q(t — t') function in the action after the expansion of the exponential. The 'external 
legs' are cut off, and this fact is indicated with the two parallel lines. 
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FIG. 3: One-loop diagrams for the two point vertex function IVj . The 'external legs' are cut off, and this fact is indicated 
with the two parallel lines. The X indicates the insertion of the composite operator <f><f>. 



for all n. The parameter U2 will also be irrelevant, and 
only the mass parameter r renormalizes. 

Consider now the renormalization group (RG) param- 
eter space spanned by r, U2, and the . The DP fixed 
point occurs at a point on the U2 axis, and the GEP fixed 
point is somewhere in the hyper-dimensional space. Since 
all the coupling constants are irrelevant in the case d > 6, 
the RG flows bring the system from the neighborhood of 
the DP fixed point towards the Gaussian DP fixed point, 
which corresponds to the DP free-field theory. Therefore 
if close to the DP fixed point we turn on a small perturba- 
tion of immunization in a system, the qualitative nature 
of the system will not change, and the epidemic process is 
expected to spread according to a mean-field DP process 
controlled by the probability of first infections. 

For 4 < d < 6, the coupling constant A^ is the most 
relevant. Consequently the higher order terms in the ex- 
pansion of the exponential in Eq. (JSJ can be neglected. 
In this case the RG flows bring the system away from 
the neighbourhood of the DP fixed point to the GEP 
fixed point which is stable. Thus, we expect the sys- 
tem to undergo an ordinary crossover from DP to the 
dynamic percolation universality class at a certain typi- 
cal time scale. However, recently it has been suggested 
that dangerous irrelevant operators may possibly lead to 
a non-trivial critical behavior different from dynamical 
percolation |2Sl |. 

To take the fluctuation effects into account, we de- 
fine the renormalized A^ as A^j and we apply standard 
methods of the perturbative renormalization group |3l| . 
In particular we make use of an e = 4 — d expansion 
around the DP critical point. The non-Mar kovian mod- 
ification in the DP action can be written then as an ad- 



ditional term of the form, 



A (1) dtd d x e(t-t')4>(x,t)^(ic,t)(f>(x,t'). (n) 



The non-locality in time, expressed by the function in 
Eq. 1111 is represented by a dashed line propagator in the 
Feynman diagrams. The renormalized coupling constant 
A^ can be determined in terms of the vertex function 
T (see Fig. 01 , evaluated at some specified normaliza- 
tion point, NP, setting a momentum scale fx. The ' in the 
vertex function's notation in Eq. 1111 is used to indicate 
that this function is calculated in the DP with immuniza- 
tion theory and is different from the one calculated in the 
ordinary DP. Near the upper critical dimension the ul- 
traviolet divergences are absorbed in the renormalization 



constant Z^, where 4>R — Z, 



and Z<h = Zi. 



A (1) - 

A R ~ 



-V ~ 

R4><i><i> 



NP 



(12) 



NP 



The structure of the corrections to the vertex function 
r 1 ! , , suggests a correspondence with the RG theory of 

DP away from the critical point, where is considered the 
renormalization of a nonzero mass term r which couples 
with the composite operator <p(x, t)(f>(x, t). One loop cor- 
rections to the corresponding two point bare vertex func- 
tion , in this theory, are depicted in Fig. [3] The nor- 
malization condition implies that the renormalized two 
point vertex function r R ^ evaluated at the normaliza- 
tion point NP, is equal to 1. 



NP 



= 1. 



(13) 



NP 
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Comparing these corrections with those appearing in the 

expression to r~ , it is possible to see that the insertion 

<t>4>4> 

of the composite operator <f><p is equivalent to the "inser- 
tion" of the dashed-continuous-line leg in Fig. [21 By in- 
spection of all possible Feynman diagrams, we conclude 
that this is valid to all orders in perturbation theory. 
Therefore the Feynman integrals involved in both A' 1 ' 
and r renormalizations are identical, and we can write: 



-A (1) rr, = r' 7j 



(14) 



Consequently from Eq. $T2&, JT3JI and (|T4*|) . it can be 

proved that, 



To describe how A^ flows under renormalization, it is 
necessary to define a Callan-Symanzik beta function. 
The dimensionless coupling constant corresponding to 



(15) 



Ag> is 



9r 



(i) 



dT 



R \ „-2-e/2 
It 

A« 



^)»- 2 -' /2 zl /2 Z^- D 2 , (16) 



where e = 4 — d, and Dp — DqZo is the renormalized 
diffusion coefficient. Then the beta function, which gives 
the differential renormalization group flow equation of gp 
is 

Mr) =^ = Sat" 2 - \ + ^ + 7* - 2 7D ], (17) 
where 7^, 7^ and are defined as 

dlnZj, 



14> = M" 



dlnZ lA 



dlnZ 



Id = M" 



(20) 



At the DP fixed point these gamma functions are related 
to the critical exponents by 7J = — <i, 7-^ = z — 
and 7^ = z — 2, and thus 

= — 08 - f|| - 1).9a + 0(g%). (21) 



The renormalization group eigenvalue y\a) j corre- 
sponding to AW, is equal to 



1 



2/ACi) 



(1 



18). 



(22) 



Since this expression is always positive, the first term 
of the expansion in Eq. 10 is relevant to all orders of 
perturbation theory around d = d c = 4. 

However, for d < d c — 4 the coupling constants are 
increasingly relevant in the expansion J!|J), and it is no 
longer meaningful to expand the exponential function. 
The system is driven away from the vicinity of the DP 
fixed point towards a non-trivial fixed point of order 4— d. 
To study this case, we apply a scaling analysis. Assume 
that position and time scale as [a;] = fc _1 and [t] = k~ z . 
In the DP theory away from criticality, the non-zero mass 
term scales as [r] — k 1 / u± . It couples to the opera- 
tor (jxj), which has a scaling dimension x^, such that, 

[4>4>] — k x i"t>. Using the fact that the action has to be 

dimensionless, we find that xi A = d + z — . Let us 

focus now on the scaling analysis of the non-Markovian 
term added to the DP action. First we consider the term 
corresponding to n = 1. AW scales as [A^] = k v ^ (1) . No- 
tice that the operators <j>(x, t)c/)(x, t) and <f>(x, t') are eval- 
uated at different space-time points, and consequently 
they are not correlated. So, each of them scales with 
its own scaling dimension, xi, and x^, — (3/v±, such 

that, [4>(x,t)(f>(x,t)] = k x io, [(f>(x,t')} = k x <f. From 
the dimensionless nature of the action, we thus obtain 
yx = —-^-{f3 — v\\ — 1). Generalizing this for higher-order 
terms in the exponential expansion we find that the scal- 
ing exponent for each \^ are, 



1 / 

— (1 + n(Mi 



/?)) 



(23) 



Since uu — (3 is always positive in d < 4 dimensions, the 
(lg) terms of the expansion JSJ are still increasingly relevant. 

Furthermore, the scaling invariance of the non- 
Markovian term can only be established if the exponen- 
tial function and its arguments are dimensionless. There- 
(19) fore the two couplings w and A scale separately with dif- 
ferent scaling exponents, i.e., [A] = k Vx and [w] = k Vu> . 
From Eq. \\1'.'>\\ we obtain, 



1 

V\ = — , 



Vu 



(24) 



Consequently the coupling constant w, is also relevant 
under renormalization. 



III. LOW FIRST INFECTION PROBABILITY 
LIMIT 



Then, as we are approaching the infrared limit (fc — ► 0), 
g' R increases, being relevant to all orders in perturbation 
theory. 



In the previous section using scaling arguments, we 
showed that the non-Markovian term is relevant for 



V 

time 



FIG. 4: Spreading in the limit of a very small first infection probability po- The figure shows a surviving cluster in 1+1 
dimensions. In this case the domain of infected/immune is hole-less and grows very slowly with time. Its boundaries can be 
considered as stationary absorbing walls. 



d < d c — 4. We can now argue that a critical spreading 
process on the horizontal line in Fig.^ is therefore driven 
away from the vicinities of the DP critical point to the 
points Aor B, where the probability of first infections is 
either very low or very high. In order to understand the 
dynamic behavior in these limits we study the asymp- 
totic spreading properties keeping the second infection 
probabilities at the critical value of DP. We consider an 
initial state with a single infected site at the origin in a 
non-infected environment. For simplicity we will focus 
on the 1 + 1-dimensional lattice. In this case there are 
not empty sites inside the epidemic clusters, which are 
composed of infected and immune sites only. 



A. Quasi-static approximation 

Let us consider the limit of a very small first infec- 
tion probability. In this regime, although first infections 
hardly ever happen, there still exists a surviving criti- 
cal epidemic process as is shown in Fig. The active 
domain is composed by re-infected sites among immune 
ones, and can be considered as bounded by almost rigid 
"walls" of empty sites. Let us define L(t) to be the dis- 
tance at time t between these two walls, such that there 
is at least one site active. If we keep the second infec- 
tion rate critical, the dynamic epidemic process can be 
thought as an effective critical DP process evolving in a 
finite system of size L(t). 

A quasi-static approximation is based on the assump- 
tion that the time scale on which L(t) grows is much 
larger than the correlation time of the DP process. 
Therefore, L(t) can be assumed constant. It is known 
that in a finite-size system with a constant width L, the 
survival probability of a DP process decays exponentially 
as P(t) ~ exp(— t/L z ), where z = Then, in a quasi- 
static finite system the survival probability can be ap- 
proximated as P s (t) ~ exp(— £/£|i), i.e. 



dP a {t) 
dt 



-aP s {t)L(t)- z . 



(25) 



where a is a non-universal amplitude factor. We need 
now to know how L(t) grows. Since the two boundaries 
are absorbing, we can apply the theory of DP with ab- 
sorbing walls . According to these results, the density 



of active sites next to the walls generated by the surviving 
clusters scales as 



Ps (t) = bL(t)-^ , 



(26) 



where b is a non-universal amplitude factor which de- 
pends on the value of po- fi s — 0.73371(2) is the sur- 
face critical exponent of DP in 1+1 dimensions. Clearly, 
dL/dt is proportional to the frequency by which the sys- 
tem attempts to infect sites at the boundary and thus 
also to the density of active sites next to the boundary: 



dL(t) 
dt 



= p oP s(t) = Po bL(t)-P>/^ 



(27) 



Therefore, the size of the domain grows as 

L(t) = (poKl+A/^)*) 1/(im/l,i) - (28) 

Inserting this result into Eq. I|25|) and solving the differ- 
ential equation we obtain 



lnP,(t) 
where a 



L{z)' z dt 



(29) 



0.947167. Hence the survival decays 
asymptotically as a stretched exponential of the form 

P s (t) = P s (0)exp(-Apo Q t 1 - Q ) . (30) 

with A — jzr^{f^) a ■ The above result implies that the 
average survival time T is finite and scales as 



T 



a/(l-a) 

■Po ■ 



(31) 



It is interesting to compare these results with recent 
findings for random walkers between movable reflectors, 
where the survival probability was shown to deca y a s 
a power law with continuously varying exponents |33| . 
We note that this result is not in contradiction with 
the present work since it would correspond to the limit 
a -► 1. 



B. Numerical simulations 

The lattice model considered is a 1+1 dimensional tri- 
angular lattice. We simulate a directed bond percolation 
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FIG. 5: The logarithm of the domain size L(i) in the problem 
of directed percolation with immunization for small values of 
po as function of time t. The numerical results are compared 
with the quasic-static approximation predictions. 



process with a first infection rate po different from a rate 
for subsequent infections or second infection rate p. 

Let us first consider the limit of a very small first in- 
fection probability, where it is difficult for the process 
to infect sites at the boundary that have never been 
infected before (see Fig. 0J. Bonds inside the infected 
region are open with the critical probability p = p c = 
0.6447001(1) while bonds leading to healthy sites 
at the boundary are open with a different probability of 
first infection po- However, as po is very small, conven- 
tional seed simulations are not suitable since most of the 
runs terminate after a very short time. For example, for 
p = 0.1 most of the runs survive for less than 100 time 
steps. Consequently we apply enrichment methods in or- 
der to circumvent this limitation. First of all, we apply 
a very simple enrichment method, in which we consider 
one lattice system. We keep the activity artificially alive 
at any time t. This is achieved by ignoring the updates 
which lead to the inactive state. Averaging over many 
realizations we measured the domain size L(t) and the 
surface density of active sites. We find that the expo- 
nents predicted by the quasi-static approximation are in 
good agreement with the values obtained with the nu- 
merical simulations as shown in Fig. [5] and Fig. The 
nonuniversal factor b in Eq. I|28|l depends on the value of 
Po- For, po — 0.01 it takes the value b = 1.65. 

To determine the survival pro bability P s (t), we apply 
another enrichment method 35], which leads to a con- 
siderable improvement. The simulation starts with an 
ensemble of N — 65 536 independent systems. Whenever 
the number of active runs becomes smaller than N/2 the 
ensemble is duplicated by creating identical copies of all 
the remaining active states. The new systems are labelled 
according to their ancestors. The survival probability P s 




i 

FIG. 6: The logarithm of the surface density of active sites 
of directed percolation with immunization for small values of 
po as function of time t. The numerical results are compared 
with the quasic-static approximation predictions. 

then is reduced by a factor of nf/rii, where rii is the 
number of initial active systems and n/ is the number 
of remaining active systems before the duplication. This 
process may be repeated as long as the ensemble has a 
sufficiently large number m of independent ancestors at 
t = 0. Using this method we were able to extend the 
temporal range of the simulation by four orders of mag- 
nitude up to t = 10 6 . Our results are shown in Fig. for 
various values of po. As can be seen, the survival prob- 
ability plotted in a double logarithmic representation is 
not a straight line, proving that it does not follow a power 
law. For the case po = 0.01, the parameter a in Eq. i)25|) 
takes the value a = 1.00(2). 



C. Lattice effects 

The scaling arguments in Section IIII Al are developed 
assuming the existence of a well defined mean value for 
the domain size L(t). To justify this assumption, we de- 
fine P(L,t) to be the probability of having a domain 
of size L at time t. Then the survival probability is 
P s (t) = J P(L,t)dL. In this subsect ion we show that 
L{t) is peaked around the mean value L(t). 

We begin by writing down a heuristic discrete master 
equation for the temporal evolution of P(L,t): 

P(L,t+l) = e- aL ~ z P{L,t)+b Po {L-l)-P°^ x 

P(L -l,t) -bpoL~^ /u ^P{L,t). (32) 

The first term in the right hand side of Eq. (|32|l describes 
the change in P{L,t) due to terminating runs according 
to Eq. H25f) integrated over the time interval [t, t+1]. The 
second term corresponds to the probability of creating 
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FIG. 7: The logarithm of the survival probability P s (t) of 
directed percolation with immunization for small values of 
po. The inset shows the corresponding number of independent 
ancestors as a function of time (see text). For po = 0.01 the 
duplication method reaches its limit since only 11 independent 
ancestors are left. 
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FIG. 8: Probability G(L, t) = f ™ P(L,t)dL computed by 
numerical integration of the master equation Eq. 13211 . The 
plots correspond to a = 0.01 and po = 0.01. The inset shows 
the scaling plot when G(L,t) scales according to Eq. 1351 . 



a domain of size L from a domain of length [L — 1). 
The third term describes a loss term corresponding to 
L — ► L + 1. We neglect the second order contributions 
+b 2 pl(L~l)- 2 ^/^P(L-2,t) and -b 2 p^(L- l)- 3 W"->- 
which only affect the behaviour at small time and small 
L regime. 

It is straightforward to show that this master equation 
is consistent with Eq. (|25|l . by summing over all L and 
replacing L by its average value L(t). On doing this 
we obtain dP s (t)/dt = —aL(t)~ z P s (t), that is the same 
equation as Eq. 1)25(1. 

Consider Eq. (|32[) with po, a and b fixed at certain 
values. Then the differential equation which corresponds 
to the master equation, Eq. (|32|) . can be written as 

dt L z 

Pob^r \L-P°/^P(L,t)] , L»l. (33) 

The solution for Eq. (jSHl is 

1> - p b((3 s /v ± + l)t] , (34) 

where L ^> 1. Consider the integrated distribution prob- 
ability G(L, t) = J™P(L,t)dL. Then G(L,t) has the 
functional form 



G(L,t) 



e -c[p b(p s /u ± +i)tr x 



L"±- - p b{ h l)t 



(35) 



with c = a/(pob(l — z + f3 s /v±)) and k = (1 — z + 
(3 s /u x )/(l + (3 8 /v x ). 

To observe the qualitative behaviour of P(L, t) and 
G(L,t), we carry out numerical integrations of Eq. I|32|l . 
As initial conditions we choose the state in which there 
is only one site infected at time t = 0. In Fig. [HI G(L, t) 
is shown for several values of time t. It behaves like a 
propagating front with an overall exponential factor of 
time t. The fact that all the G(L,t) curves intersect at 
the same point (see inset in Fig. |SJ , implies that there 
is a well defined peak and a mean value L(t) for P(L,t). 
The qualitative behaviour of G(L, t) is independent of the 
particular values of a and po. Consequently, to perform 
these calculations, we use values of a and po such that 
the effect of the exponential decay in Eq. i|35|) is not so 
pronounced. 

We now compare the predictions of the master equa- 
tion proposed in this subsection, with the results ob- 
tained from a simulation of the model. In Eq. (|32|l . we 
fixed the value of po = 0.01 and vary a. The survival 
probability calculated from direct numerical integration 
of Eq. I132[) . describes well the corresponding simulation 
result for p = 0.01, in the case a = 1 (see Fig. |HJ|. 



IV. HIGH FIRST INFECTION PROBABILITY 
LIMIT 



Let us finally consider the limit of a very high first 
infection probability po — ► 1. Again we restrict to the 
1+1-dimensional case, where the region of immune sites 
does not contain healthy sites inside. Because of the en- 
hanced spreading probability at the boundaries, the do- 
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space 



FIG. 9: Survival probability calculated using the master 
equation Eq. \YA'2}i . We set po = 0.01, a = 1 and b — 1.65. 
The result is in good agreement with the Monte Carlo simu- 
lations for po — 0.01. 



main grows rapidly by sudden avalanches of successive 
first infections (see Fig. I10[) . Since the avalanches at the 
left and the right boundary are expected to be uncorre- 
cted, it suffices to study the propagation of one of the 
boundaries. As in the previous section, we propose sim- 
ple scaling arguments in order to describe the growth of 
the infected/immune domain and the survival probabil- 
ity. 



A. Independent avalanche approximation 

In 1+1-dimensional directed bond percolation with p = 
p c and po —> 1 an avalanche is caused by a sequence of 
open bonds at the boundary. Therefore, the avalanche 
size £ is distributed exponentially as 



P(0 = (l-Po)p§ 



The quantity 

£ = - 



i 



bxpo 



1 - Pa 



Po 



1. 



(36) 



(37) 



is the average distance by which the avalanche advances 
the boundary in space. After each avalanche the pro- 
cess continues to evolve as an ordinary critical DP pro- 
cess inside the immune domain until it terminates or re- 
turns to the boundary where it releases a new avalanche. 
Thus the spreading behavior is mainly determined by the 
distribution of waiting times r between the avalanches. 
We argue that the distribution of waiting times between 
avalanches is related to the problem of local persistence 
in DP j3|| (f° r a recent review on persistence see 
e.g. [l^). The local persistence probability R(t) is de- 
fined as the probability that a randomly selected site in 




FIG. 10: Spreading in the limit of a very high first infection 
rate. 



an ordinary critical DP process starting from a homo- 
geneous initial state has not been reactivated until time 
t. It was shown that this quantity decays algebraically 
as R(t) ~ t~ e , where = 1.50(1) is the so-called lo- 
cal persistence exponent |3(| . In the present problem 
the situation is similar: Each avalanche creates locally a 
quasi-homogeneous state. The process then evolves as an 
ordinary critical DP process inside the infected/immune 
region until the boundary is revisited for the first time in 
order to release a new avalanche. However, unlike persis- 
tence studies in 1+1 dimensions, where a persistent site 
can be activated independently from the left and from 
the right, the boundary sites in the present problem can 
be infected only from one side. Hence the probability 
that the next avalanche has not yet been released decays 
as r~ e / 2 . Thus we conjecture that the waiting times 
between avalanches are distributed algebraically as 

-i-e/2 



P(r) 



(38) 



Next, we argue that these waiting times may be inter- 
preted as directed Levy nights in time j^. After each 
flight the domain size L(t) grows on average by the mean 
avalanche size £. This type of growth may be described 
by the equation 



Df /2 L{t)~l 



(39) 



0/2. 

where D t is a fractional derivative denned through its 



) 0/ 



2 c iut 



Simple 



action in Fourier space D~ t ! e wt = (iui 
dimensional analysis leads to the result 

L(t) ~ f t»/ a . (40) 

In order to compute the survival probability, we note that 
during the waiting time the outermost active site departs 
from the boundary with an average distance £(t) ~ t 1 '*. 
Obviously, the process can only terminate if this distance 
is of the same order as the size of the infected/immune 
domain L(t). Therefore, we expect the distribution of 
waiting times between avalanches to be cut off by a max- 
imal waiting time 



L z (t)~£ z t 



:0/2 



(41) 



Consequently, the probability that the process terminates 
between two avalanches is given by 



Pn = 



J?dTP(r) 



-e/2 



(42) 
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On the other hand, the cutoff due to terminating runs 
implies that the average waiting time is finite and scales 
as 



J[—drP(r) 



0/2 



(43) 



Therefore, the average loss of the survival probability per 
unit time is given by Pq/t ~ l/r max , i.e. 



dPsjt) 
dt 



-P.(f)/r a 



-P s {t)C Z t 



z ,-z@/2 



(44) 



Solving this equation, the asymptotic behavior of the 
survival probability is not a stretched exponential and is 
given by a 



P s (t) = icxp l+AC't 1 ' 



Z0/2 



(45) 



where N and A are unknown constants. In contrast to 
the previous case in Eq. 1(30(1 the exponent 1 — zG/2 ~ 
—0.185 is negative. Consequently the survival probability 
tends to a constant P s (oo) > 0, meaning that a finite 
fraction of runs survives for infinitely long time. For very 
large po this constant is expected to scale as 

l-P.(oo) oc (l- Po ) 2e/2 - (46) 



B. Numerical results 



In order to verify these results numerically, we devel- 
oped an especially optimized Monte Carlo algorithm for 
large values of pq. In order to compute the survival 
time for a run and to see how the boundaries advance 
for a given realization of open and closed bonds, it is in 
most cases not necessary to construct the entire cluster. 
Rather it suffices to construct its branches next to the 
boundaries. For example, in Fig. large parts of the 
cluster are irrelevant for the advancement of the bound- 
aries and the survival time. For this reason our algo- 
rithm first constructs the branches of the cluster next 
to the boundaries. If these branches terminate before 
the simulation time is reached, we recursively constructs 
the omitted branches in the interior of the clusters. If 
however they do not terminate, then the recursive con- 
struction is not performed. Using this technique we could 
extend the simulation time by two decades to 10 4 time 
steps, two decades less than in the previous case of low 
Po- 

Our numerical results are shown in Fig. ^] The left 
panel shows the survival probability as a function of 
time in a directed bond percolation process with po — 
0.85,0.90,0.95, and 0.98. The positive curvature of the 
lines indicates that there is no power-law scaling. Al- 
though the upward curvature is in agreement with the 
expected result <|45ll . the simulation time is not large 



enough to confirm the behaviour of Eq. Q45[) quantita- 
tively, mainly because the constant A is not known. In 
order to substantiate the assumptions made in the pre- 
vious subsection, we first measured the growing domain 
size L(t). As shown in Fig. II lb . the measured slope seems 
to tend to the predicted slope 0/2 = 0.75(1). Moreover, 
we estimated the logarithmic derivative of the survival 
probability, in which the unknown prcfactor A drops out 
(see right panel of Fig. Ill|) . Although this data set is 
quite noisy, we observe a rather clean power law. The es- 
timated slope 1.12(10) is consistent with the theoretical 
prediction £~ ze / 2 — i -1 - 185 , supporting the assumptions 
which led to the result 145(1 . 



CONCLUSIONS 



We analyzed the effects of immunization as a small per- 
turbation on the DP model and studied in detail the scal- 
ing behaviour of the theory around the DP critical point. 
We derived by an alternative method the field theoretic 
action for the model. A non-Markovian term is added 
to the DP action because of the presence of immuniza- 
tion. In the field theory, the probability for subsequent 
infections, which is different from the first infection prob- 
ability, is a function of position and time. Nevertheless, 
we assumed in the lattice model that the susceptibility 
for spreading changes only after the first infection, re- 
maining constant thereafter. This assumption does not 
changes the final results. 

The phase diagram (cf. Fig. ^) comprises two phase 
transition lines, namely, a horizontal line, where the pro- 
cess in reinfected regions shows the critical behavior of 
DP, and (in more than one spatial dimension) a curved 
transition line, where the critical behavior corresponds to 
the general epidemic process studied in Both lines 
meet at the DP critical point p = pq = p c . 

The non-Markovian term turns out to be relevant for 
d < 6. We considered the RG flows and argued that 
any point in the neighbourhood of the DP critical point 
will be driven away from it. In particular, we focused 
our study on the horizontal phase transition line with a 
critical reinfection probability p — p c . The system along 
this line, is driven away from the DP critical point as 
soon as the immunization effects are turned on as a small 
perturbation. The asymptotic behavior is determined by 
the limits of very low and high first infection probability 
(close to points A and B in Fig.^). 

We proposed simple scaling arguments for the be- 
haviour of the survival probability in both limits. We re- 
lated the corresponding exponents with those exponents 
of the critical ordinary DP theory and DP near a wall. 
The survival probability is found to obey a stretched ex- 
ponential behaviour in the low first infection probability 
limit, and it decays to a constant in the high first infec- 
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FIG. 11: Avalanche approximation in comparison with numerical results, (a) Numerically determined survival probability 
for po — 0.85, 0.9, 0.95, and 0.98 from bottom to top. (b) Numerically determined domain size L(t) for po = 0.9 (solid line) 
compared to the theoretical prediction 14011 (dotted line), (c) Logarithmic derivative of the survival probability for po = 0.9. 



tion probability limit, taking the form 

fexp(-at+ 0528 ) forpo 
s[ > K \eM+bt-° 185 ) for Po 



■ 

■ 1 



(47) 



The numerical simulations in 1 + 1 dimensions support 
these theoretical predictions, ruling out the possibility of 
asymptotic power-law scaling in both limits. The ques- 
tion as to whether the stretched exponential decay of the 
survival probability persists in higher dimensions is still 
open. 

Finally, we comment on a possible connection of these 
results with the problem of infinitely many absorbing 
states. A typical example of such systems is the pair con- 
tact process 2 A — > 3A, 2A — > 0, in which soli tary pa rti- 
cles are not allowed to diffuse In Ref. [ill IMI^ E| 
this model was studied with an effective Langcvin equa- 
tion and it was inferred that the survival probability de- 
cays as a power law at the transition point, with contin- 
uously varying ex pon ents. However, the field theoretical 
action studied in [R|, |2J, |26| is identical to the one stud- 
ied in this paper. Hence, we conclude that, according 
to the analysis presented in this paper, seed simulations 



of the pair contact process, at least in 1 + 1 dimensions, 
should show an stretched exponential behaviour and not 
a power low decay at the transition point. 
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